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... Thus logics and mathematics in the central nervous system, when viewed as lan- 
guages, must structurally be essentially different from those languages to which our 
common experience refers. ... whatever the system is, it cannot fail to differ consider- 
ably from what we consciously and explicitly consider as mathematics. 

— John Von Neumann [53] 



Abstract - Criticality in the cortex emerges from the seemingly random interaction 
of microscopic components and produces higher cognitive functions at mesoscopic and 
macroscopic scales. Random graphs and percolation theory provide natural means to de- 
scribe critical regions in the behavior of the cortex and they are proposed here as novel 
mathematical tools helping us deciphering the language of the brain. 



1.1 Introduction 

Von Neumann emphasized the need for new mathematical disciplines in order to under- 
stand and interpret the language of the brain [S3]. He has indicated the potential directions 
of such new mathematics through his work on self-reproducing cellular automata and fi- 
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nite mathematics. Due to his early death in 1957, he was unable to participate in the 
development of relevant new theories, including morphogenesis pioneered by Turing [5D] 
and summarized by Katchalsky [39] . Prigogine [56] developed this approach by modeling 
the emergence of structure in open chemical systems operating far from thermodynamic 
equilibrium. He called the patterns dissipative structures, because they emerged by the 
interactions of particles feeding on energy with local reduction in entropy. Haken [34] 
developed the field of synergetics, which is the study in lasers of pattern formation by par- 
ticles, whose interactions create order parameters by which the particles govern themselves 
in circular causality. 

Principles of self-organization and metastability have been introduced to model cogni- 
tion and brain dynamics [3B], [3D], [35] . Recently the concept of self-organized criticality 
(SOC) has captured the attention of neuroscientists [4]. There is ample of empirical ev- 
idence of cortex conforming to the self-stabilized, scale-free dynamics of the sand pile 
during the existence of quasi-stable states [T] [21 EH] ■ However, the model cannot produce 
the emergence of orderly patterns within a domain of criticality [441 132] . Bonachela [TpJ 
describe SOC as pseudo-critical and suggest to complement self-organization with more 
elaborate, adaptive approaches. 

Extending on previous studies, we propose to treat cortices as dissipative thermody- 
namic systems that by homeostasis hold themselves near a critical level of activity that 
is far from equilibrium but steady state, a pseudo-equilibrium. We utilize Freemans neu- 
roscience insights manifested in the hierarchical brain model: the K (Katchalsky) sets 
[391 [18] . Originally, K-sets have been described mathematically using ordinary differential 
equations (ODE) with distributed parameters and with stochastic differential equations 
(SDE) using additive and multiplicative noise [111113] . This approach has produced signifi- 
cant results, but certain shortcoming have been pointed out as well j 191[2~01 i22[. Calculating 
stable solutions for large matrices of nonlinear ODE and SDE that closely correspond to 
chaotic ECoG activity are prohibitively difficult and time-consuming to model on both 
digital and analog platforms [57] • In addition, there are unsolved theoretical issues in con- 
structing solid mathematics with which to bridge the difference across spatial and temporal 
scales between microscopic properties of single neurons and macroscopic properties of vast 
populations of neurons [211 123] . 

In the past decade, neuropercolation approach has proved to be an efficient tool to 
address the above shortcomings by implementing K-sets using concepts of discrete math- 
ematics and random graph theory [41]. Neuropercolation is a thermodynamics-based 
random cellular neural network model, which is closely related to cellular automata (CA), 
the field pioneered by Von Neumann who anticipated the significance of CA in the con- 
text of brain-like computing [54] . The present study is based on applying neuroperco- 
lation to systematically implement Freeman's principles of neurodynamics. Brain dy- 
namics is viewed as a sequence of intermittent phase transitions in an open system with 
synchronization-desynchronization effects demonstrating symmetry breaking demarcated 
by spatio-temporal singularities 24, 48, 261 129] . 

This work starts with the description of the basic building blocks of neurodynamics. 
Next, we develop a hierarchy on neuropercolation models with increasing complexity in 
structure and dynamics. Finally, we employ neuropercolation to describe critical behavior 
of brains and to interpret experimentally observed ECoG/EEG dynamics manifesting 
learning and higher cognitive functions. We conclude that criticality is a key aspect of the 
operation of brains and it is a basic attribute of intelligence in animals and in man-made 
devices. 
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1.2 Principles of Hierarchical Brain Models 

1.2.1 Freeman K Models: Structure and Functions 

We propose a hierarchical approach to spatio-temporal neurodynamics, based on K sets. 
Low-level K sets were introduced by Freeman in the 70s, named in the honor of Aharon 
Kachalsky, an early pioneer of neural dynamics [18]. K sets are multi-scale models, describ- 
ing increasing complexity of structure and dynamical behaviors. K sets are mesoscopic 
models, and represent an intermediate-level between microscopic neurons and macroscopic 
brain structures. The basic KO set describes the dynamics of a cortical microcolums with 
about 10 4 neurons. K-sets are topological specifications of the hierarchy of connectivity 
in neuron populations. When first introduced, K-sets have been modeled using a system 
of nonlinear ordinary differential equations (ODE) 18, 22 . K-dynamics predict the os- 
cillatory waveforms that are generated by neural populations. K-sets describe the spatial 
patterns of phase and amplitude of the oscillations. They model observable fields of neural 
activity comprising EEG, LFP, and MEG. K sets form a hierarchy for cell assemblies with 
the following components [25]: 

■ KO sets represent non-interactive collections of neurons with globally common inputs 
and outputs: excitatory in KOe sets and inhibitory in KOi sets. The KO set is the 
module for K-sets. 

■ KI sets are made of a pair of interacting KO sets, both either excitatory or inhibitory 
in positive feedback. The interaction of KOe sets gives excitatory bias; that of KOi 
sets sharpens input signals. 

■ KII sets are made of a Kle set interacting with a Kli set in negative feedback giving 
oscillations in the gamma and high beta range (20-80 Hz). Examples include the 
olfactory bulb and the prepyriform cortex. 

■ Kill sets made up of multiple interacting KII sets. Examples include the olfactory 
system and the hippocampal system. These systems can learn representations and 
do match-mismatch processing exponentially fast by exploiting chaos. 

■ KIV sets made up of interacting Kill sets are used to model multi-sensory fusion and 
navigation by the limbic system. 

■ KV sets are proposed to model the scale-free dynamics of neocortex operating on and 
above KIV sets in mammalian cognition. 

K sets are complex dynamical systems modeling the classification in various cortical 
areas, having typically hundreds or thousands of degrees of freedom. Kill sets have been 
applied to solve various classification and pattern recognition problems [221 1121 142] , In 
early applications, Kill sets exhibited extreme sensitivity to model parameters which pre- 
vented their broad use in practice [T2] . In the past decade systematic analysis has identified 
regions of robust performance and stability properties of K-sets have been derived |63II38| . 
Today, K sets are used in a wide range of applications, including detection of chemicals 
|33j . classification 12, 28], time series prediction [5] and robot navigation [371 143] . 

1.2.2 Basic Building Blocks of Neurodynamics 

The hierarchical K model-based approach is summarized in the 10 Building Blocks of 
Neurodynamics [TBI \H\ 125] : 

1. Non-zero point attractor generated by a state transition of an excitatory population 
starting from a point attractor with zero activity. This is the function of the KI set. 

2. Emergence of damped oscillation through negative feedback between excitatory and 
inhibitory neural populations. This is the feature that controls the beta-gamma 
carrier frequency range and it is achieved by KII having low feedback gain. 
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3. State transition from a point attractor to a limit cycle attractor that regulates steady 
state oscillation of a mixed E-I KII cortical population. It is achieved by KII with 
sufficiently high feedback gain. 

4. The genesis of broad-spectral, aperiodic/chaotic oscillations as background activity 
by combined negative and positive feedback among several KII populations; achieved 
by coupling KII oscillators with incommensurate frequencies. 

5. The distributed wave of chaotic dendritic activity that carries a spatial pattern of 
amplitude modulation AM in Kill. 

6. The increase in nonlinear feedback gain that is driven by input to a mixed population, 
which results in the destabilization of the background activity and leads to emergence 
of an AM pattern in Kill as the first step in perception. 

7. The embodiment of meaning in AM patterns of neural activity shaped by synaptic 
interactions that have been modified through learning in Kill layers. 

8. Attenuation of microscopic sensory-driven noise and enhancement of macroscopic 
AM patterns carrying meaning by divergent-convergent cortical projections in KIV. 

9. Gestalt formation and preafference in KIV through the convergence of external and 
internal sensory signals leading to the activation of the attractor landscapes leading 
to intentional action. 

10. Global integration of frames at the theta rates through neocortical phase transitions 
representing high-level cognitive activity in the KV model. 

Principles 1 through 7 describe the steps towards basic sensory processing, including 
pattern recognition, classification, and prediction, which is the function of Kill models. 
Principles 8 and 9 reflect the generation of basic intentionality using KIV sets. Princi- 
ple 10 expresses the route to high-level intentionality and ultimately consciousness. The 
greatest challenge in modeling cortical dynamics is posed by the requirement to meet two 
seemingly irreconcilable requirements. One is to model the specificity of neural action even 
to the level that a single neuron can be shown to have the possibility of capturing brain 
output. The other is to model the generality by which neural activity is synchronized and 
coordinated throughout the brain during intentional behavior. Various sensory cortices 
exhibit great similarity in their temporal dynamics, including the presence of spontaneous 
background activity, power-law distribution of spatial and temporal power spectral den- 
sities, repeated formation of AM spatial patterns with carrier frequencies in the beta and 
gamma ranges, and frame recurrence rates in the theta range. 

Models based on ODE and SDE have been used successfully to describe the meso- 
scopic dynamics of cortical populations for autonomous robotics [431 147] . However, such 
approaches suffered from inherent shortcomings. They falter in attempts to model the 
entire temporal and spatial range including the transitions between levels, which appear 
to take place very near to criticality. Neuropercolation and random graph theory offers 
a new approach to describe critical behavior and related phase transitions in brains. It 
is shown that criticality in the neuropil is characterized by a critical region instead of 
a singular critical point, and the trajectory of the brain as a dynamical systems crosses 
the critical region from less organized (gaseous) phase to more organized (liquid) phase 
during input induced destabilization and vise versa [3111301 127] . Neuropercolation is able 
to simulate these important results from mesoscopic ECOG/EEG recording across large 
spacial and temporal scales, as will be introduced in this essay. 

1.2.3 Motivation of Neuropercolation Approach to Neurodynamics 

We utilize the powerful tools of random graph theory (RGT) developed over the past 50 
years 16,8,9 to establish rigorous models of brain networks. Our model incorporates cel- 
lular automata and percolation theory in random graphs that are structured in accordance 
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with cortical architectures. We use the hierarchy of interactive populations in networks as 
developed in Freeman K models [TS], but replace differential equations with probability 
distributions from the observed random graph that evolve in time. The corresponding 
mathematical object is called neuropercolation |41j . Neuropercolation theory provides 
a suitable mathematical approach to describe phase transitions and critical phenomena 
in large-scale networks. It has potential advantage compared to ODEs and PDEs when 
modeling spatio-temporal transitions in the cortex, as differential equations assume some 
degree of smoothness in the described phenomena, which may not be very suitable to 
model the observed sudden changes in neurodynamics. Neural percolation model is a nat- 
ural mathematical domain for modeling collective properties of brain networks, especially 
near critical states, when the behavior of the system changes abruptly with the variation 
of some control parameters. 

Neuropercolation considers populations of cortical neurons which sustain their metastable 
state by mutual excitation, and that its stability is guaranteed by the neural refractory 
periods. Nuclear physicists have used the concept of criticality to denote the threshold 
for ignition of a sustained nuclear chain reaction, i.e., fission. The critical state of nu- 
clear chain reaction is achieved by a delicate balance between the material composition 
of the reactor and its geometrical properties. The criticality condition is expressed as the 
identity of geometrical curvature (buckling) and material curvature. Chain reactions in 
nuclear processes are designed to satisfy strong linear operational regime conditions, in 
order to assure stability of the underlying chain reaction. That usage fails to include the 
self-regulatory processes in systems with nonlinear homeostatic feedback that characterize 
cerebral cortices. 

A key question is how cortex transits between gas-like randomness and liquid-like order 
near the critical state. We have developed thermodynamic models of cortex [301127) . which 
postulates two phases of neural activity: vapor-like and liquid-like. In the vapor-like phase 
the neurons are uncoupled and maximally individuated, which is the optimal condition 
for processing microscopic sensory information at low density. In the liquid-like phase the 
neurons are strongly coupled and thereby locked into briefly stable macroscopic activity 
patterns at high density, such that every neuron transmits to and receives from all other 
neurons by virtue of small- world effects [151 159] . Local 1/f fluctuations have the form of 
PM patterns that resemble droplets in vapor. Large-scale, spatially coherent AM patterns 
emerge from and dissolve into this random background activity but only on receiving a CS. 
They do so by spontaneous symmetry breaking [31] in a phase transition that resembles 
condensation of a rain drop, in that it requires a large distribution of components, a 
source of transition energy, a singularity in the dynamics, and a connectivity that can 
sustain interaction over relatively immense correlation distances. We conclude that the 
background activity at the pseudo-equilibrium state conforms to self-organized criticality, 
that the fractal distribution of phase patterns corresponds to that of avalanches, that the 
formation of a perception from sensory input is by a phase transition from a gas-like, 
disorganized, low-density phase to a liquid-like high-density phase |30| . 



1.3 Mathematical Formulation of Neuropercolation 
1.3.1 Random Cellular Automata on a Lattice 

In large networks, such as cortex, organized dynamics emerges from the noisy and chaotic 
behavior of a large number of microscopic components. Such systems can be modeled as 
graphs, in which neurons become vertices. The activity of every vertex evolves in time de- 
pending on its own state, the states of its neighbors, and possibly some random influence. 
This leads us to the general formulation of random cellular automata. In a basic two-state 
cellular automaton, the state of any lattice point x £ Z d is either active or inactive. The 
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lattice is initialized with some (deterministic or random) configuration. The states of the 
lattice points are updated (usually synchronously) based on some (deterministic or prob- 
abilistic) rule that depends on the activations of their neighborhood. For related general 
concepts, see cellular automata such as Conway's Game of Life, cellular neural network, 
as well as thermodynamic systems like the Ising model [51 1141 1521 I62| . Neuropercolation 
models develop neurobiologically motivated generalizations of cellular automata models 
and incorporate the following major conditions: 



■ Interaction with noise: The dynamics of the interacting neural populations is inher- 
ently non-deterministic due to dendritic noise and other random effects in the nervous 
tissue and external noise acting on the population. Randomness plays a crucial role 
in neuropercolation models and there is an intimate relationship between noise and 
the system dynamics, due to the excitable nature of the neuropil. 



■ Long axon effects (rewiring): Neural populations stem ontogenetically in embryos 
from aggregates of neurons that grow axons and dendrites and form synaptic con- 
nections of steadily increasing density. At some threshold the density allows neurons 
to transmit more pulses than they receive, so that an aggregate undergoes a state 
transition from a zero point attractor to a non-zero point attractor, thereby becom- 
ing a population. In neural populations, most of the connections are short, but there 
are a relatively few long-range connections mediated by long axons. The effect of 
long-range axons are similar to small-world phenomena [61] and it is part of the 
neuropercolation model. 



■ Inhibition: An important property of neural tissue that it contains two basic types 
of interactions: excitatory and inhibitory ones. Increased activity of excitatory pop- 
ulations positively influence (excite) their neighbors, while highly active inhibitory 
neurons negatively innuence(inhibit) the neurons they interact with. Inhibition is in- 
herent in cortical tissues and it controls stability and metastability observed in brain 
behaviors |40l I38| . Inhibitory effects are are part of neuropercolation models. 



1.3.2 Update Rules 

A vertex Vi G V of a graph G(V, E) is in one of two states, s(tij), and influenced via the 
edges by d(vi) neighbors. An edge from Vi to Vj, ViVj £ E, either excites and inhibits. 
Excitatory edges project the states of neighbors and inhibitory edges project the opposite 
states of neighbors, if the neighbor's state is 1, and 1 if it is 0. Vertex's state, influenced 
by edges, is determined by the majority rule; when the most neighbors are active, the 
higher a chance for the vertex to be active, and when the most neighbors are inactive, the 
higher a chance for the vertex to be inactive. At time t = s(vi) is randomly set to or 1. 
Then, for t = 1, 2, T — 1, a majority rule is applied simultaneously over all vertices. A 
vertex Vi is influenced by a state of its neighbor Vj £ N(vi), whenever a random variable 
R(vi,t) is less than the influencing excitatory edge VjVi strength, Uji, else the vertex Vi 
is influenced by an opposite state of neighbor Vj. If the edge VjVi is inhibitory, the vertex 
Vj sends when s(vj) = 1, and 1 when s(i)j) = 0. Then, a vertex Vi gets a state of the 
most common influence, if there is such, otherwise a vertex state is randomly set to or 



1, Fig. 1.1(a) 
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Figure 1.1 Illustration of update rules in 2D lattice, (a): An example of how majority 
function works, bx edge inhibits with strength tJb,x, so it sends when s(b, t — 1) = 1. Given 
the scenario, s(x,t) is most likely. |l.l(b)| Example of 2-D torus of order 4x4. First 
row/column is connected with last row/column. Each vertex has a self-influence. Right: 
2-D torus after the basic random rewiring strategy. Two out of eighty (16 x 5) edges are 
rewired or 2.5%. 



Formula for the majority rule: 



oif E /(**,«)< ^ 



s(vi,t) = < 



lif E 



or lif £ f{v i ,t) = $& 



(1.1) 



f(vj,t) 



if u jti (s(vj,t-l)=0) >R(vi,t); else 1 

1 if Uj,i(s(vj,t-i) = l) >R(vi,t); else 



s(vj,t-l) = 



if u}j t i (s(Vj , t— 1) = 0) excites; else 1 

1 if uijAs{Vj, t— 1) = l) excites; else 



< R(vi,i) < 1 
0.5 < Uj,i(s(vj)) < 1 
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1.3.3 Two-Dimensional Lattice with Rewiring 

Lattice-like graphs are built by randomly re-connecting some edges in the graphs set on a 
regular 2-dimensional grid in and folded into tori. In the applied basic random rewiring 
strategy, n directed edges from n x 2 vertices are plugged out from a graph randomly. At 
that point, the graph has n vertices lacking an incoming edge and n vertices lacking an 
outgoing edge. To preserve the vertex degrees of the original graph, the set of plugged-out 
edges is returned back to the vertices with the missing edges in the random order. An 
edge is pointed to the vertex missing the incoming edge and is projected from the vertex 
missing the outgoing edge, Fig. |l.l(b)| Such a two-dimensional graph models a KI sets 
with homogeneous population of excitatory or inhibitory neural population. 

It will be useful to reformulate the basic subgraph on a regular grid, as shown in Fig. 
|1.2| We list all vertices along a circle and connect the corresponding vertices to obtain 
the circular representation on the top right corner of Fig. |1.2| This circle is unfolded 
into a linear band by cutting it up at one location, as shown in the middle of the figure. 
Such representation is very convenient for the design of the rewiring algorithm. Note that 
this circular representation is not completely identical to the original lattice, due to the 
different edges when folding the lattice into torus. We have conducted experiments to 
compare the two representations and obtained that the difference is minor and do not 
affect the basic conclusions. An important advantage of this band representation that it 
can be generalized to higher dimensions by simply adding suitable new edges between the 
vertices. For example, in 2D we have 4 direct neighbors and thus 4 connecting edge, in 3D 
we have 6 connecting edges, and so on. generalization to higher dimensions is not used in 
this work as the 2D model is a reasonable approximation of the layered cortex. 




0-D 1-D 2-D 3-D ... n-D 

* + + x ••• 

|V|° |V|J \V\^ 

Figure 1.2 Top-left: Example of labeled 4x4 torus in 2-dimcnsion. Top-right: 
Conversion of the 2D lattice approximately into an 2-D torus where vertices are set in circle. 
Middle: 4x4 2-D torus cut an unfolded into a linear band. Bottom section: illustration of 
the neighborhoods in higher dimensions. 
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1.3.4 Double- Layered Lattice 



We construct a double-layered graph by connecting two lattices, see Fig. 1.3(a) The top 
layer Go is excitatory, while the bottom layer Gi is inhibitory. The excitatory subgraph 
Go projects to the inhibitory subgraph Gi through edges that excite, while Gi inhibitory 
layer projects toward Go with edges that inhibit. An excitatory edge from Go influence 
the vertex of Gi with 1 when active and with when inactive. Conversely, an inhibitory 
edge from Gi influences the vertex of Go with when active and with 1 when inactive. 



1.3.5 Coupling Two Double-Layered Lattices 



By coupling two double-layered graphs, we have four subgraphs interconnected to form a 
graph G = Go U Gi U G2 U G3. Subgraph Go and Gi are coupled into one double layer 
and G2 and G3 into the other. There are 3 excitatory edges between Go and G2, Go 
and G3, Gi and G2, Gi and G3, respectively, which are responsible to couple the two 



double-layered lattices, Fig. 1.3(b) 




^ v £ Go excites G\ 
^Ch9--- V -> G x 1 if 1 else 



v £ G\ inhibits Go 
~¥~"v — > Gq 1 if else 1 




Figure 1.3 Illustration of coupled layers with randomly selected connection weights 
between them. 1.3(a) Example of a double-layer coupled system that includes an excitatory 
layer Go and an inhibitory layer Gi . Within the layers the edges are excitatory; edges from 
Go to Gi are excitatory and edges from Gi to Go are inhibitory. 1.3(b) Example of two 
coupled double-layers, all together 4 coupled layers: Go and G2 are excitatory layers, while 
GI and G3 denote inhibitory layers. Each subgraph has three edges rewired. Edges from 
Gi to Go and edges from G3 to G2 are inhibitory; the other edges are excitatory. 
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1.3.6 Statistical Characterization of Critical Dynamics of Cellular Automata 

The most fundamental parameter describing the state of a finite graph G at time t is its 
overall activation level S(G, t) defined as follows: 

|v|— 1 

S(G,t)= s (^) (1-2) 

8=0 

It is useful to introduce the normalized activation per vertex, a(G,t), which is zero if 
all sites are inactive and 1 if all sites are active. In general, a(G,t) is between and 1: 

a(G ,t ) = ^M (13) 

The average normalized activation over time interval T is given as: 

T-1 

(a) = ^_ (1.4) 

T-1 

X>(G,*)-0.5| 
(a*) = _ (1.5) 

The latter quantity can be interpreted as the average mangetization in terms of sta- 
tistical physics. Depending on their structure, connectivity, update rules, and initial 
conditions, probabilistic automata can display very complex behavior as they evolve in 
time. Their dynamics may converge to a fixed point, to limit cycles, or to chaotic at- 
tractors. In some specific models, it is rigorously proven that the system is bistable and 
exhibits critical behavior. Namely, it spends a long time in either low- or high-density 
configurations, with mostly or 1 activation values at the nodes, before a very rapid 
transition to the other state [21 [35] . In some cases, the existence of phase transitions can 
be shown mathematically. For example, the magnitude of the probabilistic component uj 



of the random majority update rule in Fig. 1.1(a) is shown to have a critical value, which 



separates the region with stable fixed point from the bistable regime [41] . 

In the lack of exact mathematical proofs, Monte Carlo computer simulations and meth- 
ods of statistical physics provide us with detailed information on the system dynamics. 
The critical state of the system is studied using the statistical moments. Based on the 
finite size scaling theory, there are statistical measures which are invariant to system size 
at the state of criticality [7]. Evaluating these measures, the critical state can be found. It 
has been shown that in addition to the system noise, which acts as a temperature measure 
in the model, there are other suitable control parameters which can drive the system to 
critical regimes, including the extent of rewiring and the inhibition strength. [58j. 

Denote this quantity ui which will be a control parameter; it can be the noise, rewiring, 
connection strength or other suitable quantity. In this work, ui describes noise effects 
and related to the system noise level e = 1 — u). Distributions of a(G,u)) and a*(G,u)) 
depend on w. At a critical value of Wj, the probability distribution functions of these 
quantities suddenly change, toi values are traditionally determined by evaluating the 4th 
order moments (kurtosis, peakedness measure) 114, [2)117, 49, 50][5l], but they can also be 
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measured using the 3rd order moments (skewness measure) or u^: 

MG) = <(a _ (a))2)2 (1.6) 

-3(G) = . <y = gog|L (i.7) 

According to finite-size scaling theory, for two structurally equivalent finite graphs G and 
G' of different sizes, 114(G) 7^ it4(G') and 143(G) 7^ 113 (G') for cj 7^ k>i. However, for uj = Wi, 
144(G) = 114 (G') and 143(G) = U;j(G'). In this study, we use both the 3rd and 4th order 
momentums to determine the critical points of the system. 



1.4 Critical Regimes of Coupled Hierarchical Lattices 
1.4.1 Dynamical Behavior of 2-D Lattices with Rewiring 

Two-dimensional lattices folded into tori have been thoroughly analyzed in [35]; the main 
results are summarized here. In the absence of rewiring, the lattice is in one of the two 
possible states: ferromagnetic and paramagnetic, depending on the noise strength u). The 
following states are observed: 

■ For oj < Wo, states and 1 are equiprobable across the graph and a distribution is 
uni-modal. This state corresponds to paramagnetic regime. 

■ For ui > cjo, one state dominates; vertex states are mostly or mostly 1, and graph's 
a distribution is bi-modal. This is the ferromagnetic condition. 

■ In the close neighborhood of the critical point ojo, the lattice dynamics is unstable. 
Immediately above luq, the vertex mostly follows its majority influences. 




Figure 1.4 Illustration of the dynamics of the 2-D torus; the case of KI architecture. 
There are two distinct dynamic regimes; top left diagram: unimodal pdf with zero average 
activation (paramagnetic regime); top middle and right: the case of bimodal pdf with a 
positive activation level(ferromagnetic regime). Lower part shows U4 = ((a — (a)) 4 )/((a — 
(a)) 2 ) 2 (left) and M3 around ujq for 2 graphs of different orders. 
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Figure 1.5 Integrated view of the relationship between the proportion of rewired edges 
and critical noise e c — 1 — too, adopted from [35]. Notations l(i), 1(2), 1(3), and 1(4) mean 
that a vertex can have 1, 2, 3, or 4 remote connections, respectively. 



For ui > ujo, the vertices cease to act individually and become a part of a group. Their 
activity level is determined by the population. This regime is defined by the condition 
that the vertices project more common activation than their receive, in the probabilistic 
sense. This transition is the first building block of neurodynamics and it corresponds to 
KI level dynamics. Fig. |1.4| illustrates critical behavior of a 2-D torus with u>o ~ 0.866; 
no rewiring. Calculations show that for the same torus having 5.00% edges rewired, the 
critical probability changes to loq ~ 0.830. Rewiring leads to a behavior of paramount 
importance for future discussions: the critical condition is manifested as a critical region, 
instead of a singular critical point. We will show various manifestations of such principle in 
the brain. In the case of rewiring discussed in this section, the emergence of critical region 
is related to the fact that same amount of rewiring can be achieved in multiple ways. For 
example, it is possible to rewire all edges with equal probability, or select some vertices 
which are more likely rewired etc. In real life such as brains, all these conditions happen 
simultaneously. This behavior is illustrated in Fig. |1.5| and it has been speculated as a key 
condition guiding the evolution of the neuropil in the infant's brain in the months following 
birth 46, 45, 2 giving rise to a robust, broad-band background activity as a necessary 
condition of the formation of complex spatio-temporal patterns during cognition. 



1.4.2 Narrow Band Oscillations in Coupled Excitatory-Inhibitory Lattices 

Inhibitory connections in coupled subgraph G = Go U Gi can generate oscillations and 
multi-modal activation states 58 . In an oscillator, there are three possible regimes de- 
marcated by two critical points loo and wi: 

■ loo is the transition point where narrow-band oscillations start. 

■ uj\ is the transition point where narrow-band oscillations terminate. 

Under appropriately selected graph parameters, a(G, ui) distribution is uni-modal when 
lo < loo, bi-modal when wo < w < wi, and quadro-modal when u) > u\, Fig. |1.6| Each 
subgraph has 2304 (5%) edges rewired within the layer and 1728 (3.75%) edges rewired 
towards the other subgraph. Just below loq, vertices may oscillate if they are excited. 
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a(G, oj) values of temporarily excited graph oscillate, but return to the basal level in the 
absence of excitation. This form of oscillation is the second building block of neuron 
inspired dynamics. When ui is further increased, i.e., uio < uj < wi, a(G,u)) exhibits 
sustained oscillations even in the absence of external perturbation. When the double- 
layer is temporarily excited, it returns to the basal oscillatory behavior. Self-sustained 
oscillation is the third building block of neuron inspired dynamics and it corresponds to 
KII hierarchy. 
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Figure 1.6 Illustration of the dynamics of coupled excitatory-inhibitory layers; the case 
of KII architecture. There are three modes of the oscillator G(5%, 3.75%) made of two 2-D 
layers Go and Gi of order 96 x 96. Top: a(G, u) of the whole graph (gray line) and of 
subgraph Go (black line) at three different cj values, which produce three different dynamic 
modes. Bottom part shows distributions of a* at 3 different u values for graphs of orders 
100 x 100 and 200 x 200; only one side of the symmetric a* distribution is shown. 



1.5 Broad-Band Chaotic Oscillations 



1.5.1 Dynamics of Two Double Arrays 



Two coupled layers with inhibitory and excitatory nodes exhibit narrow-band activ- 
ity oscillations. These coupled lattices with appropriate topology can model any fre- 
quency range with oscillations occurring between ujq and wi. For example, an oscillator 
G(5%, 2.5%) = Go U Gi has 2.5% edges rewired across the two layers and 5% edges are 



rewired within each layer. G(5%, 2.5%) has ujq » 0.840 and wi « 0.860, Fig. |1.7(b)| Simi- 
larly, oscillator G(5%,3.75%) has uio w 0.845 and ui w 0.865. Each of the coupled lattices 
in this example is of size 100 x 100. When coupling two such oscillators together, we have a 
total of 4 layers (tori) coupled into a unified oscillatory system. Clearly, the two oscillators 
coupled together may or may not agree on a common frequency. We will show that under 
specific conditions, various operational regimes may occur, including narrow band oscilla- 
tions, broad-band (chaotic) oscillations, and intermittent oscillations, and various combi- 
nations of these. One of multiple ways to couple two oscillators is to rewire excitatory edges 
between subgraphs G0-G2, G0-G3, G1-G2, and G1-G3, respectively. Two coupled oscilla- 
tors have additional behavioral regimes. When oscillators G(5%, 2.5%) and G(5%, 3.75%) 
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are coupled with 0.625% edges into a new graph G(5%, 2.5%,. 3.75%, 0.625%), uio is 
shifted to lower level, Fig. |1.8| Every parameter describing a graph influences the critical 
behavior and under appropriate conditions. If two connected oscillators with different 
frequencies cannot agree on a common mode, but together they can generate large-scale 
synchronization or chaotic background activity. 

In the graph made of two oscillators, there are four critical points ujo < ui\ < ui2 < 1^3, 
separating 5 possible major operational regimes. 

■ Case of uu < uio'- The aggregate activation distribution is uni-modal and it represents 
the paramagnetic regime in lattice dynamics analogy. 

■ cjo < w < oji : The two coupled oscillators exhibit large-scale synchronization and the 
activation distribution is bi-modal. 

■ cji < uj < UJ2'- The two coupled oscillators with different frequencies cannot agree on 
a common mode, so together they generate aperiodic background activity. 

■ u)2 < oj < UJ3: only one oscillator oscillates in a narrow band and the activation 
distribution is okta-modal. 

■ uj > U3 : none of the components exhibit narrow-band oscillation and the activation 
distribution is a hexa-modal. This corresponds to ferromagnetic state. 

The above regimes correspond to the fourth building block of neurodynamics as described 
by Kill sets. 



1.5.2 Intermittent Synchronization of Oscillations in 3 Coupled Double Arrays 

Finally, we study a system made of three coupled oscillators, each with its own in- 
herent frequency. Such a system produces broad-band oscillations with intermittent 
synchronization-desynchronization behavior. In order to analyze quantitatively the syn- 
chronization features of this complex system with over 60,000 nodes in six layers, we merge 
the channels into a single two dimensional array at the readout. These channels will be 
used to measure the synchrony among the graph components. First we find the dominant 
frequency of the entire graph using Fourier transform after ensemble averaging across the 
~I0,000 channels for the duration of the numeric experiment. We set a time window W 
for the duration of two periods determined by the dominant frequency. We calculate the 
correlation between each channel and the ensemble average over this time window W . At 
each time step, each channel is correlated with the ensemble over the window W . Finally, 
we determine the dominant frequency of each channel and for the ensemble average at time 
t and the phase shifts are calculated. The difference in phases of dominant frequencies 
measure the phase lag between the channel and the ensemble. 

We plot the phase lags across time and space to visualize the synchrony between graph 
components, see Fig. |f ,9| We can observe intermittent synchronization spatio-temporal 
patterns as the function of the uu, which represents the noise component of the evolution 



rule. For low ui values the phase pattern is unstructured, see Fig. 1.9(a) At a threshold 



us value intermittent synchronization occurs Fig. 1.9(b) Further increasing ui, we reach a 
regime with highly synchronized channels. Intermittent synchronization across large cor- 
tical areas has been observed in ECoG and EEG experiments as well. The correspondence 
between our theoretical/computational results and experimental data is encouraging. Fur- 
ther detailed studies are needed to identify the significance of the findings. 



1.5.3 Hebbian Learning Effects 

According to principles of neurodynamics, learning reinforces Hebbian cell assemblies 
which respond to known stimuli by destabilizing the broad-band chaotic dynamics an 
leading to a attractor basin through a narrow-band oscillation using gamma band career 
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Figure 1.7 Critical states obtained by coupling two double layers of oscillators with 
rewiring (case of Kill), (a): Left: the graph starts to oscillate after an impulse, but the 
oscillation decays and returns to the basal level; right: oscillation is sustained even after the 
input is removed, if wo < b: < cji. (b): Critical behavior and ljo and uj\ values for two 
different oscillators, G(5%,2.5%) (up) and G(5%,3.75% (bottom). 



frequency. To test this model, we implemented learning in the neuropercolation model We 
use Hebbian learning, i.e., the weight from node i to node j is incrementally increased if 
these two nodes have the same state (1 — 1 or — 0). The weight fromi toj'j incrementally 
decreases if the two nodes have different activations (0 — 1 or 1 — 0). Learning has been 
maintained during the 20 step periods when input was introduced. 

In the basal mode without learning, the 3 double layers can generate broad-band chaotic 
oscillations, see Fif 1.10(a) where the lower row is zoomed-in version of top row. The 
spikes in Fig 1.10(a) indicate the presence of input signals. Inputs have been implemented 
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Figure 1.8 Coupling two KII oscillators into Kill, lowers ujq, the point above 
which large-scale synchronous oscillations are formed. The top plots show values for 
G(5%, 2.5%, 3.75%, 0.625%), G(5%,2.5%, 3.75%, 1.25%), and G(5%, 2.5%, 3.75%, 2.5%). 
Lower part shows typical activation of graph for uj below and above ojq. 



by flipping a few % of the input layer nodes to state 1 (active) for the duration of 20 
iteration steps. During the learning stage, inputs are introduced 40 times (20 steps each), 
at regular intervals of 500 iteration steps. Without learning, the activity returns to the 



low-level chaotic state soon after the input ceases. Fig 1.10(b) shows that a narrow- 
band oscillation becomes prominent during learning, when a specific input is presented. 
After learning, the oscillatory behavior of the lattice dynamics is more prominent, even 
without input, but the learnt input elicits much more significant oscillations. This is the 
manifestation of the 6th and 7th principles of Freemans neurodynamics, and it can be 
used to implement classification and control tasks using the neuropercolation model. 

Further steps, toward the 8th and 9th principles involve the analysis of the amplitude 
modulation (AM) pattern, which is converted to a feature vector of dimension n. Here 
n is the number of nodes in our lattice which, after some course-graining, correspond 
to the EEG/ECoG electrodes. They provide our observation window into the operation 
of the cortex or cortex simulation. The feature vector is used to describe or label the 
cortical state at any moment. The AM pattern is formed by a low-dimensional, quasi- limit 
cycle attractor after synaptic changes with learning. The synaptic weights are described 
in a matrix, and the combination of different modalities to form Gestalts is done by 
concatenation of feature vectors. Research into this direction is in progress. 



1.6 Conclusions 

We employed neuropercolation model to demonstrate basic principles of neurodynamics. 
Neuropercolation is a family of stochastic models based on the mathematical theory of 
probabilistic cellular automata on lattices and random graphs and motivated by structural 
and dynamical properties of neural populations. The existence of phase transitions has 
been demonstrated in probabilistic cellular automata and percolation models. Neuroper- 
colation extends the concept of phase transitions to large interactive populations of nerve 
cells near criticality. Criticality in the neuropil is characterized by a critical region instead 
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Figure 1.9 Dcsynchronization to synchronization transitions; (a) at uj below critical 
value there is no synchronization; (b) when uj exceeds a critical value, the graph exhibits 
intermittent large-scale synchronization. 



of a singular critical point. The trajectory of the brain as a dynamical systems crosses the 
critical region from less organized (gaseous) phase to more organized (liquid) phase during 
input induced destabilization and vise versa. Neuropercolation simulates these important 
results from mesoscopic ECOG/EEG recording across large spacial and temporal scales. 

We showed that neuropercolation is able to naturally model input-induced destabiliza- 
tion of the cortex and consequent phase transitions due to sensory input and learning 
effects. Broad-band activity with scale-free power spectra is observed for unlearnt condi- 
tions. On the other hand, Hebbian learning results in the onset of narrow-band oscilla- 
tions indicating the selection of specific learnt inputs. These observations demonstrate the 
operation of the first 7 building blocks of neurodynamics. Our results indicate that the 
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Figure 1.10 Learning effect in coupled oscillators, (a) Activity levels without learning; 
first row contains extended time series plots, 2nd row has zoomed in sections. Input spike is 
shown at every 500 steps. The activity returns to base level after the input ceases. |1.10(b)| 
Activity levels with learning;Thc oscillations are prominent during learning and maintained 
even after the input step is removed (decayed). 



introduced Hebbian learning effect can be used to identify and classify inputs. Considering 
that our model has 6 coupled lattices each with up to 10,000 nodes, massive computational 
resources have been involved in the simulations. Parallel chip implementation may be a 
natural way to expand the research in the future, especially when extending the work for 
multi-sensory Gestalt formation in the cortex. 
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